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ABSTRACT 

We present SKIRT (Stellar Kinematics Including Radiative Transfer), a new Monte Carlo 
radiative transfer code that allows the calculation of the observed stellar kinematics of a dusty 
galaxy. The code incorporates the effects of both absorption and scattering by interstellar 
dust grains, and calculates the Doppler shift of the emerging radiation exactly by taking into 
account the velocities of the emitting stars and the individual scattering dust grains. The code 
supports arbitrary distributions of dust through a cellular approach, whereby the integration 
through the dust is optimized by means of a novel efficient trilinear interpolation technique. 

We apply our modelling technique to calculate the observed kinematics of realistic mod- 
els for dusty disc galaxies. We find that the effects of dust on the mean projected velocity 
and projected velocity dispersion are severe for edge-on galaxies. For galaxies which deviate 
more than a few degrees from exactly edge-on, the effects are already strongly reduced. As 
a consequence, dust attenuation cannot serve as a possible way to reconcile the discrepancy 
between the observed shallow slopes of the inner rotation curves of LSB galaxies and the pre- 
dictions of CDM cosmological models. For face-on galaxies, the velocity dispersion increases 
with increasing dust mass due to scattering, but the effects are limited, even for extended dust 
distributions. Finally, we show that serious errors can be made when the individual velocities 
of the dust grains are neglected in the calculations. 

Key words: dust, extinction - galaxies: kinematics and dynamics - galaxies: spiral - radiative 
transfer 



1 INTRODUCTION 

During the last decade of the past century there has been a vivid 
discussion about the opacity of spiral galaxies. This discussion 
was initiated by Disney, Davies & Phillipps (1989) and Valen- 
tijn (1990), who countered the conventional view that spiral galax- 
ies are optically thin over their entire optical discs (Holmberg 1958; 
de Vaucouleurs, de Vaucouleurs & Corwin 1976; Sandage & Tam- 
mann 1981). This issue has not yet been settled, although an overall 
consensus seems to be emerging that spiral galaxies are, in gen- 
eral, at least moderately optically thick within their optical disc, 
i.e. that they have central face-on optical depth of at least unity. 
We will not repeat a detailed overview of the different points of 
view and the large variety of observational tests, and refer to Davies 
& Burstein (1995), Xilouris et al. (1997), Kuchinski et al. (1998) 
and Calzetti (2001) for excellent reviews. We Just wish to point 
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out an important observational constraint which, in our opinion, 
is not always given enough attention: the far-infrared (FIR) emis- 
sion of spiral galaxies. This emission originates from the absorption 
of starlight by dust grains, which is thermally re-emitted at longer 
wavelengths. The most recent results show that about 30 per cent 
of the bolometric luminosity of late-type galaxies is FIR emission 
by interstellar dust (Popescu & Tuffs 2002). Notice that this figure 
refers to normal quiescent spiral galaxies; special classes of galax- 
ies such as starburst galaxies emit even higher fractions of their 
energy at FIR wavelengths. These facts demonstrate quite convinc- 
ingly that interstellar dust forms an important constituent of spiral 
galaxies. 

The physical processes of absorption and scattering by dust 
grains should hence be taken into account when interpreting and 
modelling the observed optical properties of galaxies. We would 
like to stress that all optical observables will to some extent be af- 
fected by dust, including the observed kinematics. Indeed, when 
interstellar dust grains absorb or scatter photons, the kinematical 
information contained within these photons, in the form of the 
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Doppler shift, will also be absorbed or scattered. It is important 
to realize that most of the kinematical data of galaxies are optical 
measurements. Firstly, stellar kinematics are always measured from 
optical absorption lines. In principle, the attenuation effects of in- 
terstellar dust could be minimized by using the sensitive ^^CO ab- 
sorption feature at 2.29 /xm (Gaffney, Lester & Doppmann 1995). 
Measuring stellar kinematics at near-infrared (NIR) wavelengths 
is, however, observationally extremely difficult due to the high sky 
surface brightness. Secondly, the vast majority of spatially resolved 
rotation curves of spiral galaxies have been measured using optical 
emission lines (usually the Ha line). They can also be measured 
using the Hi or CO lines, where dust extinction is not an issue. 
However, these methods are observationally much more expensive 
than the more straightforward Ha observations. 

Currently, the only way to investigate the effects of dust atten- 
uation on the observed kinematics of galaxies is by means of de- 
tailed radiative transfer modelling. Many different approaches ex- 
ist to handle the radiative transfer problem, and various teams have 
adopted different techniques to study the effect of dust on the pho- 
tometry and spectral energy distributions (SEDs) of galaxies (e.g. 
Bruzual, Magris & Calvet 1988; Witt, Thronson & Capuano 1992; 
Byun, Freeman & Kylafis 1994; Wise & Silva 1996; Bianchi, Fer- 
rara & Giovanardi 1996; Corradi, Beckman & Simonneau 1996; 
Baes & Dejonghe 2001a). The effect of dust attenuation on the ob- 
served kinematics, however, is largely unexplored, and limited to 
the rotation curves of spiral galaxies. Davies (1990) argued that the 
effects of dust absorption on the observed optical rotation curve of 
a disc galaxy can lead to a severe underestimation of the true rota- 
tional velocity in the inner regions. This argument was quantified 
by Bosma et al. (1992), who calculated the effect of dust absorption 
on the observed rotation curve of edge-on galaxies. The most de- 
tailed work on this subject comes from Matthews & Wood (2001), 
who calculate the effects of dust attenuation on the rotation curve 
of low surface brightness (LSB) galaxies through Monte Carlo sim- 
ulations. 

We have embarked on a program to systematically investigate 
the effects of dust on the observed kinematics of galaxies. In pre- 
vious papers we focused on the effects of dust attenuation on the 
observed stellar kinematics of elliptical galaxies. Initially, we only 
considered absorption by dust grains in our models (Baes & De- 
jonghe 2000; Baes, Dejonghe & De Rijcke 2000). The inclusion of 
absorption in the calculation of the observed kinematics is a fairly 
simple procedure: it basically adds a weight to the contribution of 
each individual star along a line of sight. The radiative transfer 
problem becomes much more complicated however when scatter- 
ing is also included. In this case, photons can leave their original 
path, such that any star can contribute to the observed kinematics 
in any line of sight. Moreover, not only the stellar velocities but 
also the individual dust grain velocities should be taken into ac- 
count in the calculation of the observed kinematics. We argue that 
the Monte Carlo method is the only radiative transfer modelling 
method in which kinematical information can be included in an el- 
egant and straightforward way. Using a one-dimensional spherical 
Monte Carlo code, we found that the effects of dust scattering are 
fairly important: dust attenuation can serve as an additional or al- 
ternative explanation for the stellar kinematical evidence of a dark 
matter halo around ellipticals (Baes & Dejonghe 2001b, 2002a). 

In this paper, we tackle the more complicated problem of in- 
vestigating the effect of dust attenuation on the observed kinematics 
of disc galaxies. We have written a new Monte Carlo code (SKIRT, 
acronym for Stellar Kinematics Including Radiative Transfer) that 
can handle any geometry of stars and dust. This code is presented 



in section|2| and tested in section|3| We present our model and the 
results of our Monte Carlo simulations in section |4| These results 
are discussed in section|5| and section|S|presents our conclusions. 



2 DESCRIPTION OF THE SKIRT CODE 

The basic characteristics of Monte Carlo radiative transfer have 
been explained at length by various authors (e.g. Cashwell & Ev- 
erett 1959; Mattila 1970; Witt 1977; Fischer, Henning & Yorke 
1994; Bianchi et al. 1996). In essence, a Monte Carlo radiative 
transfer code follows the life of a very large number of individual 
photons. A photon is, at each stage in its existence, characterized by 
various quantities, such as its position r, propagation direction k, 
wavelength A and Doppler shift u as measured in a frame of rest.^ 
After being emitted, photons propagate on straight lines through 
the interstellar medium until they either interact with a dust grain or 
leave the galaxy. The various interactions alter the properties of the 
photon, according to random numbers generated from the appro- 
priate probability functions. When at last the photon escapes from 
the galaxy, its final properties are recorded. After recording a large 
numbers of photons in this way, the global observed properties of 
the system can be calculated. 

We will not repeat an in-depth description of the principles 
and the numerous equations of Monte Carlo radiative transfer here, 
as they are well described in the articles listed above. Instead, we 
will focus on a number of aspects which make our code different 
from the existing ones. These differences are mainly our choice of 
the dust grid and the possibility to include kinematical information 
into the radiative transfer calculations. 

2.1 The emission process 

When no kinematical information is included in the radiative trans- 
fer calculations, photons are characterized at each moment by a 
wavelength A, a position r and a direction k. Initial values for 
these quantities must be generated randomly from the emitting stel- 
lar system, which can be composed of several stellar components.^ 
The initial emission direction fco can be sampled from the unit 
sphere, the initial position ro is sampled from the spatial distri- 
bution of the stellar component, and the initial wavelength Ao is 
sampled from the SED of the stellar component at the position ro. 
The SKIRT code contains a library with a set of common stellar 
components (including spherical Jaffe, Hernquist, Plummer and de 
Vaucouleurs models, and axisymmetric exponential, sech, isother- 
mal disc models), from which it is straightforward to sample a ran- 
dom position (see Appendix A). Another library contains a set of 
SEDs, including Planck functions and realistic tabulated stellar and 
galaxy SEDs (Kinney et al. 1996; Pickles 1998). This library also 
contains so-called monochromatic SEDs, which are described by a 
Dirac delta function. Sampling a random wavelength from the latter 
is trivial, because all photons are emitted with the same wavelength. 

^ Throughout this paper, we will make no distinction between the mea- 
sured Doppler shift AA and the corresponding line-of-sight velocity u, 
which are related by the expression AA/Aq = u/c. 

^ We describe the code as adopted for radiative ti'ansfer calculations in a 
dusty galaxy. The code can, however, equally well be applied to any other 
environment. Terms like "stellar component" must therefore be interpreted 
in a broad way, and can also refer to a single star, emitting gas, an AGN, 
etc., i.e. any source of photons of high enough energy. 
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When the SKIRT code is used to calculate the observed kine- 
matics of a dusty galaxy, photons must also carry with them the 
kinematic signature of the star that has emitted them, in the form 
of a Doppler shift. Because this is a one-dimensional quantity, not 
all kinematical information of the star can be transmitted to the 
observer, and it is important to consider how the observed Doppler 
shift relates to the velocity vectors of the star and the scattering dust 
grains. This problem was discussed in detail in section 2.3 of Baes 
& Dejonghe (2002a), where it is shown that the general expression 
for the Doppler shift after M scattering events is 
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where v* and Vi^ are the velocity vectors of the emitting star and 
the ith dust grain respectively. The initial value of the Doppler shift 
attached to the photon is hence the component of the star in the 
direction of the original emission, uq = ■ ko. 

To complete the initialization of a photon, we also need to 
generate a line-of-sight velocity from the appropriate probability 
distribution. In this case, this function is the spatial line-of-sight 
velocity distribution (LOSVD) (?i*(ro, fco, u), which describes the 
probability for a star at a position ro to have a velocity component 
u in the direction feo. In general, the spatial LOSVD is a marginal 
probability distribution of the stellar distribution function F* (r , i;), 
which describes the probability density of the stellar component in 
six-dimensional phase space, 

(t)*{r,k,u)^ — 11 F*(r,w)duxi duxa, (2) 
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where n, (r ) is the stellar number density and the integration covers 
the entire velocity space perpendicular to the direction fe. Thus, in 
order to be of use for a kinematical Monte Carlo run, each stellar 
component must allow the generation of a velocity from its spatial 
LOSVDs for arbitrary r and fc. 

The spatial LOSVDs can be calculated analytically for a num- 
ber of self-consistent spherical models. Apart from isotropic mod- 
els (e.g. Plummer 1911; Henon 1959; Jaffe 1983; Hernquist 1990), 
this is possible for the anisotropic families of Plummer and Hern- 
quist models described by Dejonghe (1987) and Baes & De- 
jonghe (2002b). For more complicated dynamical models, how- 
ever, an analytical evaluation of the spatial LOSVDs is not possible. 
For such models, we approximate the velocity distribution at each 
position in the galaxy as a local gaussian distribution. This means 
that at each position r, the distribution function can be written as 
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Here (y\,V2-,V3) and (ai, (72, 0-3) represent the mean velocities 
and velocity dispersions in the directions (ei, 62, 63), the princi- 
ple axes of the velocity ellipsoid at the position r. The orientation 
of the velocity ellipsoid and the values of the mean velocities and 
velocity dispersions can vary from position to position, such that 
the distribution function ^3) represents a fairly general distribution 
function. But a distribution function of this form has the advantage 
that the spatial LOSVDs can be calculated exactly for any direction 
k by applying the formula j2j. We find that the spatial LOSVD will 
also be a gaussian distribution. 
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with parameters 

u = k\v\ + ^2^2 + ki,vz, (5) 

a'i = k\a\^k\a\ + k\a\, (6) 

where obviously (fci, ^2,^3) are the components of the vector fe 
with respect to the orthonormal reference system (61,62,63). If 
hence, for each stellar component, we specify the orientation and 
parameters of the velocity ellipsoid at each position r, we can eas- 
ily sample line-of-sight velocities from the spatial LOSVDs in an 
arbitrary direction, and hence initialize the photons' initial Doppler 
shifts. 



2.2 The dust iteration 

The life of a single photon can be thought of as a loop, whereby 
at each iteration (representing a physical process) we must update 
its position r, propagation direction fc, wavelength A and Doppler 
shift u. The initial values are determined randomly at the emis- 
sion phase, and change at every scattering event. As we only con- 
sider coherent scattering, and we are not taking into account the 
re-emission of photons at longer wavelength^, the only wavelength 
change of a photon along its path is due to the varying Doppler 
shift. These wavelength variations are so tiny that the optical prop- 
erties of the dust (extinction curve, dust albedo and asymmetry pa- 
rameter) do not vary. The wavelength of a photon can hence be 
considered fixed throughout its lifetime. We must therefore only 
update the three quantities r, fc and u at each iteration, from their 
old values (r^-i . . . ) to the new ones (r^ . . .). This proceeds in 
several steps. 

The first step in the iteration consists of determining whether 
an interaction with a dust grain will take place or whether the pho- 
ton will leave the galaxy. To do this, we sample an optical depth 
T\ from an exponential distribution and compare it to Tpath.A, the 
total optical depth along the path. When t\ > rpath,A, the photon 
will leave the galaxy; in the other case, the photon will interact 
with a dust grain. This interaction can either be a scattering or an 
absorption event, which is determined by the scattering albedo. In 
the former case, the photon will continue its journey through the 
galaxy, in the latter case, the loop is ended. 

The second step in the loop, if the photon is scattered, is the 
determination of the position of the scattering. Therefore we have 
to translate the sampled optical depth t\ to a physical path length 
s. With this path length known, the new position is Vi — Vi^i + 
sfci-i. 

The third step in the iteration is the determination of the new 
propagation direction ki of the photon. It is found by sampling a di- 
rection from the probability density p(fci) = $\{ki-i, ki), which 
represents the scattering phase function. Various phase functions 
are built into the SKIRT code, including the isotropic phase func- 
tion and the anisotropic Henyey-Greenstein phase function. The 
sampling of a direction from these phase functions can be per- 
formed analytically. 

The final step is to update the photon's Doppler shift. The rel- 
ative orientation of the propagation directions before and after the 
scattering event cause a change in Doppler shift [from equation 0] 



Ui — Ui-i = Vi. ■ [ki — ki-i) 
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3 For the present study, we are primaiily interested in the attenuation of 
starlight by interstellar dust at optical wavelengths, where the contribution 
of re-emitted photons is negligible. 
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To calculate the updated Doppler shift, we need to sample a dust 
grain velocity from the dust velocity distribution. Therefore, the 
dust components in the SKIRT library must not only be specified by 
a spatial distribution, but also by the entire phase space distribution 
Fiir, v). Analogous to the stellar distribution function, we assume 
that the dust velocity field can be described by a trivariate gaussian 
distribution. Because the dispersions in the interstellar medium are 
fairly small compared to the rotational velocities however, we re- 
strict ourselves to assuming an isotropic dispersion tensor, such that 



look for the cell in which the interaction will take place. Since the 
opacity increases linearly with path length within a single cell, this 
inversion of the equation is then straightforward. 

In the SKIRT code, a UDD grid is included based on a three- 
dimensional cartesian grid. The grid cells can be chosen arbitrarily, 
in order to achieve an efficient grid for a variety of dust distribu- 
tions. An algorithm is included to construct clumpy dust distribu- 
tions from an underlying smooth dust distribution, as described in 
Witt & Gordon (1996) and Bianchi et al. (2000). 
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with V and a a being the mean velocity vector and the velocity dis- 
persion of the dust. Similarly as in the stellar case, we don't have 
to sample a full three-dimensional velocity vector for the dust, but 
only a component in the direction 

All ki— I 



(9) 



This means that we just have to sample a velocity from a one- 
dimensional gaussian distribution with mean velocity u — v ■ k'i 
and dispersion CTu — a^. 

2.3 Integration through the dust 

The first step in the iteration described in the previous subsection 
requires the calculation of the total optical depth along a given path. 
It is found through the integral 
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where Mi\{r) represents the total (absorption plus scattering) opac- 
ity of the dust at a position r. To execute the second step in the it- 
eration, we need to convert an optical depth tx into a physical path 
length s{r,k,T\) along a given path, which is done by solving the 
equation 



K.x{r + s'k) ds' = Tx 



(11) 



for s = s{r,k,Tx). Only for a few simple geometries, such as con- 
stant density or some spherical dust distributions, can these integra- 
tions be performed analytically. In general, the calculation of the 
two quantities Tp.^th,x{i', k) and s(r, k, tx) must be performed nu- 
merically. These integrations are usually the most time-consuming 
part of Monte Carlo codes. Therefore, it is important to think about 
an efficient and flexible way to do this; we have considered two 
approaches. 

2.3.1 The UDD grid 

In order to allow a completely arbitrary distribution of dust (in- 
cluding clumpy dust distributions), most modern radiative transfer 
Monte Carlo codes adopt an approach that can be described as a 
uniform dust density (UDD) grid. It consists of dividing space into 
a number of cells, and attaching a uniform dust density to each cell, 
for example the value at the centre of the cell. This uniform den- 
sity (and hence opacity) makes it very easy to calculate the optical 
depth [equation jlOU along a given path: simply calculate the dis- 
tance that the photon runs through each cell it crosses, multiply this 
distance with the opacity in the cell, and add all the pieces together. 
The inversion of the equation <1 1> is not much more difficult. First, 



2.3.2 The IDD grid 

In addition to the UDD grid, we explore a new approach to han- 
dle the integration through the dust: an interpolated dust density 
(IDD) grid. The basis is the same: we divide space into a number 
of cartesian grid cells, but, instead of attaching a uniform opacity 
to each cell, we use the correct values of the opacity at the eight 
border points of the cells, and we apply a simple trilinear interpola- 
tion routine to determine the opacity in each of the points within a 
grid cell. Consider a photon at the position r = {x,y,z),a. position 
within the fc)'th cell with border points Xi ^ x ^ a^i+i, etc. 
We define the dimensionless quantities p, q and r as 

p= , q= , r= . (12) 

Xi+i — Xi yj+1 — yj Zk+i — Zk 

Since p increases linearly from to I as x moves linearly from Xi 
to Xi+\, and similarly for q and r, the trilinear interpolation rule for 
the opacity simply reads 

>^x(r) = (1 - p) (1 - q) (1 - r) mj^k 
+ (1 -P) (1 - q)rK,ij^k+i 
H 

+ pqr Ki+ij+i^k+i, (13) 

where Hij,k, etc. represent the opacity in the border points. Now 
consider the fraction of the path that remains in this cell. This path 
can be parametrized using the expression r + sk, with the path 
length s as a free parameter. In each of the points along this path 
(as long as we remain in the same cell), we can calculate the opacity 
by a similar formula as <I3> . 

Kxir + sk) = [1 - p{s)] [1 - q{s)] [1 - r{s)] Kij^k 
+ [1 - p{s)] [1 - q{s)] r(s) Ki,j,fc+i 
H 

+ p{s)q{s)r{s)K,i+i,j+i,k+i, (14) 

where the parameters p, q and r are now linear functions of the path 
length s, 

X S kx Xi 



p{s) 



r(s) 



Xi+l - Xi 

y + sky-yj 

Vj+i - Vj 
z + skz — Zk 



(15) 
(16) 
(17) 



Zk+l — Zk 

Along a path within a single cell, the opacity is thus a cubic poly- 
nomial in s, the distance traveled along the path, 

3 



K.x{r + sfc) = ^ 



(18) 



with constant coefficients a™.. With this expression for the opacity, 
the calculation of the optical depth Tp-ati.xij' ,k;) is very straight- 
forward. Indeed, it is sufficient to calculate, for each cell the path 
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Figure 1. A comparison of the accuracy of the integrations through a UDD 
dust grid (constant density) and the IDD dust grid (trilinear interpolation). 
The three panels contain the relative errors of the opacity k\ (r), the optical 
depth along a path Tpath,A(^i k) and the path length s(r, fc, t\) for 1000 
random photons. The dust grids in both cases contain the same cells. The 
different colors coirespond to dust grids with 100^ (green), 200^ (blue) and 
300^ (red) cells respectively. 

passes through, the two boundary positions. This enables us to cal- 
culate the path length s within the cell and the coefficients am- 
Given these coefficients, the portion of the optical depth in that cell 
is found by substitution into the quartic polynomial 



(r + s'fe) ds' = 



(19) 



Also, the calculation of the path length s{r,k,Tx) is straightfor- 
ward. First, we determine the cell in which the optical depths at the 
border points bracket the value of , in a similar way as described 
above to calculate the total optical depth along the path. The de- 
termination of s then comes down to finding the root of a quartic 
equation. The obvious way of finding this root is by means of a 
Newton-Raphson iteration, because we can easily compute both the 
function values and the derivatives of this function. This routine has 
a very fast quadratic convergence (the number of significant digits 
approximately doubles with each iteration). Usually, no more than 
two iterations are necessary. 

2.3.3 Comparison of the grids 

To compare the accuracy and efficiency of the two approaches, we 
applied both techniques to a simple spherical dust geometry, char- 
acterized by an opacity 

/ ^2X— /2 

Kx{r)(x[l + — ] . (20) 



The integrations through the dust, i.e. the calculation of the two 
functions TpMh,\{r,k) and s{r,k,T\), can be performed exactly 
for such a dust component with a = 2, 3 or 5 [see Baes (2001) for 
details]. These exact results were compared with the results from 
the UDD and IDD grids. In figureQwe plot the relative error of the 
calculated values of the opacity, the total optical depth along a path 
and the path length corresponding to a given optical depth for 1000 
randomly generated photons for such a model. For a given number 
of dust grid cells, the accuracy of the integrations in the IDD grid is 
much better than those in the UDD grid. However, in terms of com- 
putational efficiency, comparing the two grids with the same num- 



bers of cells is not fair, because the integrations in the UDD grid 
are less complicated and hence faster. We found that the IDD grid 
integrations are a factor 3 slower in the mean. On the other hand, 
in order to reach a similar accuracy as the IDD grid, many more 
cells are necessary in the UDD grid. Apart from slowing down the 
integrations, this will also increase the memory requirements sig- 
nificantly. After several tests, we found that the IDD grid is com- 
putationally preferable above the UDD grid, in particular when the 
system harbours a large dynamical range of dust densities (such 
as an exponential profile which is appropriate in disc galaxies). It 
should be noted however, that the UDD grid still remains a very 
useful alternative, in particular to include clumpy dust distributions 
into the SKIRT code. 



2.4 Optimization of the code 

The procedure described at the beginning of the previous section 
is the most basic version of a Monte Carlo iteration step. When 
executed as such, however, it would be very inefficient. Through- 
out the years a number of 'intelligent tricks' have been presented 
which increase the computational efficiency considerably. Most of 
them attach a weight to each photon, which can be altered during 
its lifetime. We incorporate three such techniques into our SKIRT 
code. 

The first optimization to the basic Monte Carlo scenario is to 
turn all interactions into scattering events. In reality, of course, the 
interaction between a photon and a dust grain can either be an ab- 
sorption or a scattering event; the nature of this event is determined 
by the albedo cox of the dust grains. The problem with absorptions 
is that the photon is lost, and therefore does not contribute to the 
observed radiation field anymore. To overcome this problem, we 
force each interaction to be a scattering event and alter the weight 
of the photon after each scattering by a factor ujx in order to com- 
pensate for the fraction of absorbed photons. 

The second technique is the peeling-off procedure introduced 
by Yusef-Zadeh et al. (1984). In the basic Monte Carlo iteration, 
each photon leaves the galaxy at a certain stage, and it can be clas- 
sified in bins according to its position and direction. This approach 
has a number of disadvantages. Firstly, if we are interested in only 
one single observing position, most of the photons are not used. 
Secondly, if the observer has the misfortunate of being located in 
a direction in which not many photons leave the galaxy, he might 
have to wait for a very long time to detect enough photons to ob- 
tain good statistics. Thirdly, if one wants to compare the view of a 
system from two observing positions close to each other, the direc- 
tion bins must be chosen to be very small, which requires a very 
large number of photons. The peeling-off technique can be used to 
tackle these problems. It consists of creating a new photon during 
each scattering [resp. emission] process, which is scattered [resp. 
emitted] exactly in the direction of the observer. This photon will be 
detected by the observer. Its weight is altered by a factor $(fc, feobs) 
[resp. 1] to compensate for the probability of scattering [resp. emis- 
sion] in the direction feobs, and a factor exp[— rpath,A(f', feobs)] to 
compensate for extinction along the line-of-sight. 

The third optimization we apply is the principle of forced first 
scattering (Witt 1977), which is useful in the relatively low opti- 
cal depth regime. A problem in such systems is that most of the 
photons leave the galaxy directly without a single interaction with 
a dust grain, and therefore many photons need to be generated in 
order to have a good statistics of the scattered radiation. To over- 
come this problem, we split the emitted photon into two parts. The 
first part of the photon directly leaves the galaxy and obtains a 
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weight exp(— Tpath.A) (because we adopt the peehng-off technique, 
we don't have to consider this fraction anymore). The other part of 
the photon with the remaining weight is forced to be scattered at 
least once, which is achieved by sampling an optical depth tx such 
that T\ < Tpaih.A- This technique is called the principle of forced 
first scattering, but nothing prevents us applying it only to freshly 
emitted photons: it can equally well be applied to scattered photons 
(Bianchi et al. 1996). The number of forced scatterings we apply 
in our calculations is usually 3, but this value can be set arbitrarily 
according to the studied system. 

2.5 The detection and data reduction processes 

The last phase in the Monte Carlo cycle is the detection phase, 
whereby the information contained in the photon is transmitted to 
the observer Because of the peeling-off procedure, the number and 
position of the observers can be chosen randomly. At each observ- 
ing position, the user can choose from several simulated instru- 
ments to detect the photons: photometers, low-resolution and high- 
resolution spectrographs. Each instrument consists of a number of 
pixels on the plane of the sky and they come in two geometries: 
rectangular and circular. Each pixel contains a detector, the nature 
of which depends on the kind of instrument. 

The photometers are the most simple instruments, whereby 
the detector in each pixel is just a counter. Every time a photon 
hits the pixel and its wavelength is within the detector's bandwidth, 
its weight is added to the intensity. At the end of the simulation, 
the 2D surface brightness image of the system is directly obtained. 
The photometers are designed such that each photon remembers 
the number of scatterings it has undergone, and this information is 
also used by the photometric detectors. This allows us to construct 
images for each individual scattering (direct emission, photons that 
have been scattered once, etc.), which can be useful to disentangle 
the effects of absorption and scattering. 

The low-resolution spectrographs work in a very similar way, 
except that there is not just one detector, but a range of detec- 
tors corresponding to different wavelength bins. In this way a 
data cube is produced (actually a set of data cubes, one for ev- 
ery scattering), very similar to the data cubes obtained by radio or 
Fabry-Perot observations. The high-resolution spectrographs basi- 
cally produce similar data cubes, with Doppler shift bins replacing 
the wavelength bins. Combining the different Doppler shift bins at 
a given position x on the sky, we obtain the LOSVD or line profile 
4>p{x, u), which describes the entire projected velocity distribution. 
From the LOSVDs, the 2D mean projected velocity field Vp{x) and 
projected velocity dispersion field a"p(a;) can be calculated. 

From a conceptual point of view, high-resolution and low- 
resolution spectrographs work in a very similar way. Both separate 
the incoming photons in different bins according to the extra infor- 
mation carried by the photon, which is the wavelength in the for- 
mer and the Doppler shift in the latter case. From a computational 
point of view, however, there is an important difference between the 
two. Indeed, the wavelength of the photon is actively used through- 
out the Monte Carlo iteration, because the optical properties of the 
dust grains depend on the photon's wavelength. On the contrary, a 
photon's path does not depend on its Doppler shift. Moreover, the 
changes in Doppler shift experienced by a photon when it moves 
through the dust are independent of the initial stellar Doppler shift. 
Therefore, we can optimize our code in the following way. Instead 
of initializing the photon's Doppler shift by generating a uq from 
the stellar spatial LOSVD, we postpone this contribution of the 
stellar velocity until later and initialize uo to zero. When the pho- 



ton leaves the galaxy after Af scattering events, it will then carry a 
Doppler shift 



Um = ^ ^d; • (fei — fci-l), 



(21) 



where we still have to add the stellar Doppler shift in order to ob- 
tain the correct value. Instead of generating a single initial stel- 
lar Doppler shift uo from the spatial LOSVD 0, (ro, fco, it), we 
split the photon over all velocity bins. Each fraction will then be 
weighted by a weight factor that represents the probability that 
the final Doppler shift of the photon will fall within the corre- 
sponding bin. In order to fall within a bin with boundaries Uk and 
Uk+i, the photon must have an initial stellar Doppler shift satisfy- 
ing Uk ^ Uo + iiji/ ^ Uk+i. The weight factor corresponding to 
the fcth bin will hence equal 



Wk 



(j),{ro,ko,u) du. 



(22) 



When we substitute the gaussian spatial LOSVD ^4} into this ex- 
pression, we can evaluate this weight factor exactly 



Wk = - 



erf 



y-k+i 



erf 



(23) 



This method is of the same kind as the optimization techniques de- 
scribed in section l2!4l and also increases the efficiency of the SKIRT 
code in a significant way. 



2.6 General architecture and performance 

The SKIRT code is implemented in the ANSI standard C-n- lan- 
guage and makes optimal use of the object-oriented structure of 
this language. The output of the code consists of FITS files, which 
are created by means of the FITSIO and CCFITS packages. 

The run time of the SKIRT code with reliable output depends 
critically on the number of photons, the number of observing di- 
rections, the spatial and velocity resolution of the instruments and 
the number of cells in the dust grid. For the simulations described 
in section|4] we found that about 5 x 10^ photons are necessary 
to obtain accurate results. The code spends 60 per cent of the time 
moving through the dust grid, while most of the remaining time 
is spent in the calculation of the spatial LOSVDs. When run on a 
modern Pentium IV PC, the run time of a single simulation (with 7 
observing positions) is of the order of 50 to 100 hours. 



3 TESTING THE CODE 

The code was first run without dust included in order to test the pro- 
jection algorithm and the calculation of the observed kinematics. 
We used simple non-rotating spherical geometries which allow an 
analytical expression of the surface brightness and projected veloc- 
ity dispersion profiles (Dejonghe 1987; Hemquist 1990; Jaffe 1983, 
Baes & Dejonghe 2002b). Also the surface brightness profiles of 
face-on and edge-on exponential discs (which can be calculated 
analytically) are accurately recovered. Moreover, when we attach 
a simple velocity structure to an exponential disc, with a flat ro- 
tation curve and exponentially decreasing velocity dispersions, the 
observed kinematics can be also calculated analytically in the face- 
on and edge-on directions. These are recovered very well, proving 
the correct behaviour of the code. 
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and do not include a bulge, because we are primarily interested in 
the disc kinematics. 




Figure 2. The B-band surface brightness distribution of an edge-on disc 
galaxy, for various values of the optical depth. The model is similar to the 
BT0.3 model described by Byun et al. (1994), i.e. an exponential stellar 
disc with scale length /i* = 4 kpc and scale height = 0.35 kpc, a 
de Vaucouleurs bulge with = 1.6 kpc and a bulge-to-total luminosity 
ratio of 0.3. The dust is also exponential with scale length hj = 4 kpc and 
2(j = 0.14 kpc, and the face-on central optical depths in the different plots 
ai'e rv = 0, 0.5, 1, 2, 5 and 10 (from top to bottom). Surface brightness 
contours are shown with a step of = 1 mag/kpc^. 



Next, we compared the results of the SKIRT code with our 
previous Monte Carlo calculations of the observed kinematics of 
elliptical galaxies (Baes & Dejonghe 2001b, 2002a). The Monte 
Carlo code used in these papers also allowed for the calculation 
of the observed kinematics, but only for spherical models, chosen 
such that the integrations through the dust could be performed an- 
alytically. We found an excellent agreement between these results 
and those calculated with the new SKIRT code. 

Finally, we compared our (photometric) results with the radia- 
tive transfer calculations of other authors. In figure |2| we plot the 
two-dimensional edge-on surface brightness distribution of an Sb 
galaxy model in the B band for various values of the optical depth. 
The models are the same models as the BT0.3 models of Byun et 
al. (1994), and they consist of a stellar disc, a stellar bulge and a 
dust disc. Comparing this figure with their figure 4b, we find an 
excellent agreement between the two results. Notice that Byun et 
al. (1994) adopted a completely different approach to solve the ra- 
diative transfer problem. 



4 THE OBSERVED KINEMATICS OF DUSTY DISC 
GALAXIES 

We apply the SKIRT code to investigate the effect of dust absorp- 
tion and scattering on the photometry and the stellar kinematics of 
normal galactic discs. We limit our modelling only to a stellar disc 



4.1 Presentation of the model 

Our model consists of a simple axisymmetric double exponential 
disc galaxy, with stellar emissivity 



■ exp 



hi 



exp 



(24) 



The luminosity Lx is a free scaling parameter in our models, as 
both the input emissivity and the output light profile scale with 
the luminosity. For the scale length and scale height we chose the 
values adopted by Byun et al. (1994) to describe the Milky Way, 
/i, = 4 kpc and = 350 pc. As we focus on the observed kine- 
matics, which are always measured in a very narrow wavelength 
region, we chose a trivial monochromatic SED centered on the V 
band. 

In order to run a kinematical simulation for such a disc galaxy, 
we also need to provide details on the kinematical structure of the 
disc. We assume that the stellar velocities at each position in the 
disc can be represented by a trivariate Gaussian distribution. 



F-t{r, v) oc exp 



2a| 



2a2 



2a1 



(25) 



Such a velocity distribution was first proposed by Schwarzschild 
(1907) to describe the distribution of stellar velocities in the solar 
neighbourhood. The mean azimuthal velocity and the three velocity 
dispersion components can vary from position to position. The en- 
tire kinematical structure of the disc will be completely determined 
if we have expressions for the four functions v^p(r\ aR{r), (Ttp{r) 
and (Jzir). Since the system is axisymmetric, these quantities will 
depend only on the cylindrical radius R and the height z. More- 
over, because stellar discs are thin and the contribution of photons 
high above or under the symmetry plane of the galaxy will be quite 
small due to the vertical exponential decrease, we will assume that 
these four functions are a function of R only. 

For the profile, we assume the arctan profile, a simple two- 
parameter fitting function that fits the observed rotation curves of 
disc galaxies fairly well (e.g. Courteau 1997), 



2w„ 



■ arctan 



(26) 



We have adopted an amplitude v^y^^ — 220 km/s and a scale radius 
Ro = 0.5 kpc, which results in a rotation curve that has a steep ini- 
tial gradient and reaches its flat part fairly quickly. The internal ve- 
locity dispersion profiles in disc galaxies are not well constrained. 
For a thin self-gravitating disc, the Jeans equations yield that the 
vertical velocity dispersion scales with the square root of the sur- 
face density. If we assume a constant mass-to-light ratio, we obtain 
the following parametrization. 



(r j = azo exp 



R 
2/i* 



(27) 



As all dispersions in a stellar disc are thought to exist from the 
gradual stellar heating of an initial cold disc, it makes sense to 
assume that the velocity dispersions in the other direction are in- 
timately coupled to the vertical dispersion. Therefore, we assume 
a constant axis ratio of the velocity ellipsoid in our model, such 
that gr and a,p have the same spatial variation as Gz- Notice that 
such behaviour is found to be in agreement with the sparse exist- 
ing data on the velocity dispersions in external galaxies (Gerssen 
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et al. 1997, 2000). The only remaining quantities that need to be 
set are normalizations of the dispersions. We adopted the values 
((T_Ro, o^o, azo) = (100, 75, 50) km/s. At a galactocentric radius 
of 8 kpc, this gives velocity dispersions of (37,28,18) km/s, which 
roughly corresponds to the measured velocity dispersions in the so- 
lar neighbourhood. 

The dust was chosen to be distributed in a similar way as the 
stellar distribution. 



..(r) = — expf-- 



exp 

Zi 



(28) 



whereby t\ (here) represents the face-on optical depth through the 
centre of the galaxy. 



K,x{r)dLZ. 



(29) 



For the dust scale length and scale height, we again choose the val- 
ues from Byun et al. (1994), /ij = 4 kpc and Zi = 140 pc. Notice 
that the dust is hence significantly more confined to the plane of 
the galaxy, a feature which gives rise to the well-known dust lanes 
in edge-on galaxies. For the velocity field of the dust grains, we 
assumed an isotropic gaussian distribution. 



Fi{r, v) oc exp 



? + vl 



2al 



(30) 



We assumed the same mean velocity as the stellar distribution, and 
a velocity dispersion = 8 km/s, which is a typical value for 
the dispersion of the cold gas in the interstellar medium (Sofue & 
Rubin 2001). Finally, for the optical properties of the dust we took 
the average Galactic values (Gordon, Calzetti & Witt 1997). 

These functions and parameters completely determine our 
model, except for the one free parameter rv, the central face-on op- 
tical depth in the V band (the corresponding optical depth at other 
wavelengths is determined by the assumed optical properties of the 
dust). As discussed in the Introduction, the optical depth in spiral 
galaxies is still a matter of debate. We ran our models for the op- 
tical depth values rv = 0, 0.5, 1, 2, 5 and 10, in order to cover a 
wide range of possible scenarios. 

4.2 Modelling results 

In this section we present and describe the results of our Monte 
Carlo simulations. We will not discuss the effects of dust on the 
photometry, since such results have been discussed at length by 
other authors, such as Byun et al. (1994), Corradi et al. (1996) and 
Bianchi et al. (1996). Instead we concentrate on the effects of dust 
on the observed kinematics, in particular on the mean projected 
velocity Vf,(x) and the projected velocity dispersion Of,{x). 

In figures|3|and|4]we compare, at various inclination angles, 
the 2D mean projected velocity and projected velocity dispersion 
fields of an optically thin galaxy (left) with the corresponding fields 
of a dusty galaxy with optical depth rv = 1 (centre) and rv = 5 
(right) respectively. Notice that the mean projected velocity field 
is not plotted in the face-on direction, as there is no rotation per- 
pendicular to the galactic plane. To facilitate the interpretation of 
these 2D kinematical fields, we plotted in figure |5| the major and 
minor axis kinematical profiles of our models at various inclination 
angles, explicitly as a function of the optical depth. 



4.2.1 Edge-on galaxies 

For edge-on galaxies, the effects of dust on both the mean projected 
velocity and the projected velocity dispersion are very strong, as 



Table 1. The effect of dust attenuation on the slope of the mean pro- 
jected velocity curve. This slope is expressed through the major axis ra- 
dius at which the rotation curve reaches half of the asymptotic value 
Dp™" = 220 sin i kms~^. These radii ai'e tabulated for the various incli- 
nations and optical depths considered in our simulations. For a infinitely 
thin and dust-free galaxy, the intrinsic value would be 0.5 kpc, independent 
of the inclination angle. 



rv 


i = 90 


i = 88 


i = 85 


i = 80 


i = 60 





2.05 


1.58 


1.24 


0.96 


0.66 


0.5 


4.25 


2.22 


1.35 


0.95 


0.66 


1 


5.69 


2.93 


1.69 


1.04 


0.65 


2 


7.26 


3.71 


2.10 


1.24 


0.65 


5 


9.23 


4.77 


2.72 


1.56 


0.70 


10 


10.63 


5.58 


3.21 


1.82 


0.75 



is obvious from the bottom row panels of figure |3| and |4] and the 
panels in the right column of figure |5| The projected mean ve- 
locity profile along the major axis rises more and more slowly as 
the optical depth increases, and tends towards a solid body rota- 
tion. These trends were expected, and in agreement with the ar- 
guments of Davies (1990) and the absorption-only calculations of 
Bosma et al. (1992). To quantify this change of slope, we have 
calculated, for each optical depth, the major axis radius at which 
the mean projected velocity reaches half of its maximum value 
^max _ sin i km s~^ . These values are tabulated in the second 
column of tableQ and demonstrate that the effects of dust attenu- 
ation are very severe, even for small optical depths. The reason for 
this shallower slope of the rotation curve is that we mainly see the 
radiation from the optically thin part of the galaxy, which is more 
and more restricted to the outer part of the galaxy as the optical 
depth increases. 

In the same way, we can understand the effects of dust at- 
tenuation on the projected velocity dispersion profile, where also a 
prominent dust lane is visible in the contour plot: we mainly see the 
light (and hence the Doppler shift) of the stars from the outer part 
of the disc, where the intrinsic velocity dispersion is much smaller. 

4.2.2 Galaxies at intermediate inclinations 

Bosma et al. (1992) found that the signature of dust absorption on 
the rotation curve depends strongly on the inclination angle. For 
inclinations which deviate as few as 5 degrees from purely edge-on, 
they found that the effects of dust on the rotation curve is already 
severely reduced, even for optical depths up to rv ~ 10. This result 
remains the same if scattering is also taken into account, because 
the absorption effects of dust strongly dominate in highly inclined 
galaxies. TableQquantitatively demonstrates these results. For the 
projected velocity dispersion, a similar effect is valid, at least for 
the major axis dispersion profile. When the entire projected velocity 
dispersion field (or just the minor axis dispersion profile) is taken 
into account, the asymmetrical signature of the dust can still be 
recognized for inclinations down to 80 degrees. 

Moving to intermediate inclinations, the effects of dust attenu- 
ation on the observed kinematics quickly become very limited. For 
example, there is hardly any effect noticeable for an inclination of 
60 degrees, which can easily be seen by comparing the correspond- 
ing panels in figure|3|and|4| 
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Figure 3. The effect of dust attenuation on tlie projected mean velocity field of our galaxy disc models. Each single panel is a diagram with contours of equal 
projected mean velocity (the so-called spider diagram). The total field of view along the major axis is 40 kpc, i.e. 5 scale lengths at each side of the galaxy 
centre. The different contours in the panels are drawn with a step Avp = 20 km/s. The five rows correspond to different inclination angles, ranging from 
i = 60 (top row) over i = 80, 85, 88 to the edge-on i = 90 (bottom row). The three different columns correspond to different values of the optical depth: the 
left column is the dust-free model, whereas the middle and right ones con'espond to ry = 1 and tv = 5 respectively. 



4.2.3 Face-on galaxies 

In the light of the previous results, we would expect that the effects 
of dust attenuation on the observed kinematics of face-on galaxies 
are completely neglegible. And indeed, on the mean projected ve- 
locity there is no effect at all: the projected mean velocity is zero 
for face-on galaxies, regardless of whether dust is present or not. 
However, we do see a conspicuous effect of dust attenuation on the 
projected velocity dispersion: it increases with increasing optical 
depth. This cannot be an absorption effect, because all stars along a 
given line of sight in our face-on galaxy discs have the same veloc- 
ity dispersion in the direction of the observer. The effect arises 
as a result of photons emitted in the central regions in a direction 
parallel to the plane of the galaxy. When these photons are scat- 
tered into a vertical direction, they have a large probability to leave 
the galaxy. They will then contribute the high-velocity information 
of the stars in the central regions to the LOSVDs, causing high- 
velocity wings and hence an increased velocity dispersion. 

The effect is similar to what we found for the projected veloc- 
ity dispersion profile of elliptical galaxies containing a diffuse dust 
component (Baes & Dejonghe 2002a). Compared to these models, 
where the signature of the dust is very significant, it is fairly modest 
for the face-on discs. The reason for this weaker impact is the fact 
that the velocity dispersion profile in our disc galaxy models have 
a fairly shallow slope, compared to the elliptical galaxy models of 
Baes & Dejonghe (2002a). In disc galaxies, the Doppler shifts con- 
tributed by the scattered photons will in the mean differ less from 
the Doppler shifts contributed by the direct photons as in elliptical 
galaxies. 

We could ask ourselves the question whether this conclusion 
is biased by the chosen star-dust geometry of our disc galaxy mod- 
els. Indeed, Baes & Dejonghe (2002a) found their results depended 
strongly on the relative star-dust geometry: the scattering effect was 
strong only for models in which the dust was more extended than 
the stars. In our disc galaxy models, the dust distribution has the 
same scale length as the stars, whereas the scale height is smaller. 



such that the dust is always sandwiched between the stars. Recently 
however, evidence has been found for extended distributions of 
cold dust in spiral galaxies, both from the radiative transfer mod- 
elling of the dust lanes in edge-on galaxies (Xilouris et al. 1997, 
1998, 1999) and from ISO and SCUBA dust emission (Alton et 
al. 1999; Davies et al. 1999). In order to investigate whether this 
would bias our results, we considered a new set of models with 
an extended dust distribution. They were constructed by increas- 
ing the dust scale length or the dust scale height by factors up to 
three, while keeping the face-on optical depth fixed. Running new 
face-on simulations for these models, we found that the strength of 
the scattering effect on the projected velocity dispersion increases 
only minimally, with deviations of the order of a few per cent. The 
modest effects found in figure |5| are hence representative for disc 
galaxies. 



5 DISCUSSION 

5.1 Stellar kinematics in disc galaxies 

One of the most useful and promising applications of stellar kine- 
matics in disc galaxies is their contribution to mass modelling. In 
its simplest form, this mass modelling consists of decomposing the 
observed rotation curve into contributions by the disc and the dark 
halo (assuming the contributions of a bulge and molecular/atomic 
gas are negligible). Even in this most simple scenario, there is still 
a virtually complete degeneracy: the mass-to-light ratio of the disc. 
Authors have often adopted a solution in which the contribution 
of the disc is maximized, leading to the so-called maximum disc 
hypothesis (van Albada & Sancisi 1986). There are however, no 
convincing arguments to favour a particular value for the stellar 
mass-to-light ratio, and often scenarios in which the dark matter 
halo dominates the disc are equally possible (e.g. Lake & Feinswog 
1989). Stellar kinematics can, in principle, provide the necessary 
additional constraint to break this degeneracy. Indeed, for a self- 



10 Baes et al. 




Figure 4. The effect of dust attenuation on the projected velocity dispersion tield of our galaxy disc models. Each single panel is a diagram with contours of 
equal projected velocity dispersion. The different contours are drawn with a step Aup = 5 km/s. The layout of this figure is similar to figure|3] except that an 
extra panel at the top is added representing the face-on direction. 



gravitating disc, the vertical velocity dispersion scales linearly with 
the square root of the surface density of the disc, and can therefore 
be used to constrain the mass-to-light ratio of the disc. 

The most obvious candidates for such an analysis would be 
face-on galaxies, because their projected velocity dispersion di- 
rectly reflects the vertical velocity dispersion in the disc. We found 
that for a dusty disc galaxy, the observed velocity dispersion is 
affected by scattering. The effects are modest, however, even for 
models with an extended dust distribution. Alas, face-on galaxies 
do have strong observational disadvantages. Foremost, their low 
surface brightness make it very hard to obtain the stellar kinemat- 
ics out to sufficiently large distances to isolate the disc component. 
Moreover, the shape and amplitude of the rotation curve cannot be 
measured directly. The amplitude can be estimated fairly accurate 
using the TuUy-Fisher relation, whereas the shape could be approx- 
imated through the (controversial) universal rotation curve formal- 
ism (Persic, Salucci & Stel 1996). This is unsatisfactory, however, 
because there is a large spread in universal rotation curves, and the 
details of mass modelling usually depend crucially on small scale 
features of the rotation curve which statistical methods cannot in- 
corporate. 

Edge-on galaxies don't suffer from these observational prob- 



lems: their surface brightness is much higher, and their apparent 
rotation curve can be determined directly. However, the interpreta- 
tion of the kinematical data is non-trivial. Firstly, the vertical ve- 
locity dispersion must be linked to the observed dispersion, which 
is a mixture of contributions from the radial and tangential velocity 
components. Based on the observed axis ratio of velocity ellipsoid 
in the solar neighborhood, Bottema (1993) estimated the vertical 
velocity dispersion in edge-on galaxies by means of the radial ve- 
locity dispersion evaluated at one scale length. It has been shown, 
however, that the velocity ellipsoid axis ratios in disc galaxies can 
vary substantially (Gerssen et al. 2000). Secondly, our modelling 
indicates that both the projected mean velocity and the projected 
velocity dispersion field in edge-on galaxies are strongly affected 
by interstellar dust. 

Galaxies of an intermediate inclination (z < 80°), however, 
are most suitable for such an analysis. Gerssen et al. (1997) have 
shown that for galaxies with an intermediate inclination, the three 
components of the velocity ellipsoid can, in principle, be deter- 
mined by studying the kinematics along major and minor axes. 
Also, because the rotation curve of inclined galaxies can directly 
be determined and their surface brightnesses are high enough to 
allow the measurements of their kinematics in a reasonable time. 
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Figure 5. The observed kinematics of our galaxy disc models along the major and minor axes. The panels on the two rows represent the mean projected 
velocity profiles and projected velocity dispersion profiles along the major axis, and the bottom panel show the minor axis projected velocity dispersion 
profiles. The different columns correspond to different inclination angles, ranging from face-on (left) to edge-on (right). Notice that the projected velocity 
profiles are the observed profiles, and are not corrected for inclination. The different curves in each panel correspond to different values of the optical depth: 
TV = (black), TV = 0.5 (green), ry = 1 (magenta), ry = 2 (cyan), ry = 5 (blue) and -ry = 10 (red). 



they form ideal targets for such a study. Moreover, determining the 
three components of the velocity ellipsoid is very useful to con- 
strain the dynamical history of disc galaxies, because the various 
proposed mechanisms leave a different imprint on the velocity el- 
lipsoid. Our modelling demonstrates that the observed kinematics 
of inclined disc galaxies can be directly interpreted and are not bi- 
ased by dust attenuation. We argue that studying the observed kine- 
matics of a substantial set of optically smooth inclined galaxies 
could thus greatly contribute to understanding the mass structure 
and dynamical history of spiral galaxies. 

5.2 Optical rotation curves in disc galaxies 

The results of our modelling are important concerning the deriva- 
tion of the rotation curve of spiral galaxies. Whereas the outer parts 
of the rotation curves of spirals were the preferred field of interest 
for nearly three decades, it is the inner slope of the rotation curve 
which nowadays stands in the spotlight. To measure the slope in 
this inner region accurately, a high spatial resolution is required, 
which makes the usual 21 -cm measurements less suitable, except 
for the most nearby galaxies. A useful alternative are CO rotation 
curves, which can achieve a high resolution in both the spatial and 
velocity directions, but these observations are quite costly in ob- 



serving time. The main alternative remains optical Ha observa- 
tions, either long-slit or Fabry-Perot. 

Using Ha rotation curves, several authors (Blais-Ouellette, 
Amram & Carignan 2001; Borriello & Salucci 2001; de Blok, 
McGaugh & Rubin 2001; de Blok & Bosma 2002) have recently 
found shallow slopes for the inner rotation curve of a significant 
number of (mainly LSB and dwarf) galaxies. Such slopes are not 
in agreement with the results of high-resolution cosmological N- 
body simulations, where it is found that dark matter halos have a 
strong cusp, and hence that rotation curves of dark matter domi- 
nated galaxies should be steeply rising (Navarro, Frenk & White 
1997; Moore et al. 1997). A number of "solutions" have been pro- 
posed to explain this discrepancy, amongst them the suggestion that 
the internal extinction in the centre of galaxies could be a major 
source of uncertainty. Our models show that, indeed, an edge-on 
galaxy with an observed slowly rising rotation curve could in prin- 
ciple be a dusty galaxy with an intrinsically steep rotation curve. 
We found, however, that this effect cannot be obtained with realis- 
tic amounts of dust if the galaxies deviate more than a few degrees 
from exactly edge-on, confirming the results of Bosma et al. (1992) 
and Matthews & Wood (2001). In addition, because LSB galaxies 
are assumed to be dust poor (McGaugh 1994, Tully et al. 1998; 
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Matthews et al. 1999), we have demonstrated that dust effects can- 
not explain the observed discrepancies. 

We do however demonstrate that dust effects are important for 
galaxies within a few degrees from edge-on, and that dust attenua- 
tion must be seriously taken into account in such cases. 

5.3 The importance of the dust grain velocities 

The calculation of the final Doppler shift of a photon in the SKIRT 
code takes into account the velocity information of the star and the 
individual velocity vectors of each individual scattering dust grain. 
This means that for every scattering event, we need to sample a 
velocity from the dust velocity field, and this is a costly operator. 
It would be very convenient from a computational point of view, 
if we could neglect the dust grain velocities. We have done this in 
our previous Monte Carlo simulations (Baes & Dejonghe 2001b, 
2002a), where we studied the observed stellar kinematics of ellip- 
tical galaxies. In these models, both stars and dust grains are sup- 
ported by random motions, and because the dust is colder than the 
stars, the extra Doppler shifts due to scattering are generally smaller 
than the stellar Doppler shift. In a spiral galaxy, however, this ar- 
gument does not hold anymore, because the velocities of both stars 
and dust grains are dominated by the rotation and are therefore of 
the same magnitude. 

In order to investigate the importance of including the indi- 
vidual dust grain velocities into the code, we ran the same mod- 
els again but without taking the dust velocities into account. In 
figure |6| we plot the major and minor axes mean projected veloc- 
ity and projected velocity dispersion profiles of our models calcu- 
lated this way. Comparing this figure with the corresponding this 
figure with figure |5| we find that the effects of dust on the mean 
projected velocity profiles are significantly overestimated when the 
dust grain velocities are not taken into account, in particular for in- 
clined galaxies. The differences between the projected velocity dis- 
persion curves in both figures are even worse: when the dust grain 
velocities are neglected, the apparent effect of dust attenuation is 
to increase the dispersion at nearly all radii and all inclinations, by 
factors of up to 200 per cent or more for realistic optical depths. 

To understand this behaviour, consider a photon emitted by a 
star with velocity in the direction feo, and assume this photon 
is scattered by a dust grain with velocity va into the direction fei, 
whereafter it manages to escape the galaxy. When we do not take 
the dust grain velocity into account, the observed Doppler shift of 
the photon will simply be ito = ■ fco. On the contrary, the correct 
observed Doppler shift of the photon reads 

tti = ■ fco + Vi ■ (fei - feo) (31) 
= {v, - Vd) ■ ko + Vd ■ ki. (32) 

In a disc galaxy, the motion of both stars and dust grains is domi- 
nated by the same rotation. If the pathlength of the photon between 
its emission and scattering is short, the velocity vectors of the star 
and dust grain will therefore be very similar, i;, ~ Vi, such that 
Ml « i;* • fci. The new correct Doppler shift is hence a typical line- 
of-sight velocity of a star in the new propagation direction fci of 
the photon. Therefore, taking the dust motion into account more or 
less forces the dust grain to adopt a Doppler shift that is "appropri- 
ate" for the new propagation direction. Along a given line-of-sight, 
the velocity information carried by scattered photons will, in the 
mean, only modestly deviate from the velocity information carried 
by photons directly emitted in same direction. If we don't take the 
dust grain velocities into account, the photon can contribute, to the 



observed kinematics in a given direction, a line-of-sight velocity 
typical for a completely different direction. 

In the elliptical galaxy models from Baes & Dejonghe 
(2002a), this does not make much difference, because the velocities 
of the stars have similar magnitudes in the different directions. In a 
rotating disc galaxy however, this makes a huge difference because 
the stars have hugely different velocities in different directions: the 
vertical and radial motions are determined by the dispersion and are 
of the order of a few tens of km/s, whereas the azimuthal motion 
is dominated by the rotation and is of the order of 200 km/s. If the 
dust grain velocities are not correctly taken into account, scattered 
photons can thus cause erroneous extremely-high-velocity wings 
in the LOSVDs. These give rise to modest but significant errors in 
the mean projected velocity, and very large errors in the projected 
velocity dispersion. 

As a remark, we want to point out that Matthews & Wood 
(2001) found that "also including the Doppler shifts arising from 
the relative bulk motion of the scattering dust particles had a neg- 
ligible effect on the final rotation curves". Whereas this statement 
is indeed true for the mean projected velocity of highly inclined 
galaxies, we nevertheless judge that this could easily be misinter- 
preted as a justification to neglect dust velocities altogether. It must 
be realized that for the calculation of the entire LOSVDs, and for 
the projected velocity dispersion in particular, the dust velocities do 
play a crucial role, and must be taken into account. 



6 CONCLUSIONS 

In this paper we have presented a novel Monte Carlo code that can 
take velocity information into account, and can hence be used for 
kinematical modelling of dusty objects. We applied this code to 
calculate the observed kinematics of dusty disc galaxies. The main 
results of this paper can be summarized as follows. 

(i) A correct inclusion of kinematical information into radiative 
transfer problems requires the inclusion of the velocities of both the 
emitting stars and the scattering dust grains. We see no other way 
of tackling this problem except with Monte Carlo techniques. 

(ii) A new approach is presented to optimize the integration 
through the dust in a Monte Carlo code, using a trilinear interpo- 
lation instead of a constant opacity within each cell. We compared 
both kinds of grids, and find that, to obtain a similar accuracy, the 
new approach is more efficient in terms of computation time and 
memory. 

(iii) The effects of dust attenuation on the kinematics of edge-on 
disc galaxies are severe. Both the mean projected velocity and the 
projected velocity dispersion are severely affected, even for modest 
optical depths. Therefore we strongly advise to always interprete 
the stellar kinematics and optical rotation curves of edge-on galax- 
ies very cautiously. 

(iv) For galaxies which are more than a few degrees from strictly 
edge-on, the effects of interstellar dust on the observed kinematics 
are much weaker. Therefore, we argue that dust attenuation cannot 
be invoked as a possible mechanism to reconcile the observed slope 
of LSB galaxies with the predicted CDM cosmological models. 

(v) Dust attenuation does not affect the kinematics of intermedi- 
ately inclined galaxies. Such galaxies hence form the ideal targets 
for stellar kinematical studies, in particular to constrain the mass 
structure and to study the kinematical history of disc galaxies. 

(vi) The projected velocity dispersion of face-on galaxies in- 
creases slightly due to scattering of dust grains into the line of sight. 
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Figure 6. The major and minor axes observed kinematics of our disc galaxy models, calculated without taking the dust velocities into account. The layout is 
similar to figurelsl 



The effects are relatively small, however, even for extended dust 
distributions. 

(vii) Neglecting the extra Doppler effect caused by the scatter- 
ing medium results in incorrect projected kinematics, in particu- 
lar if the velocity components of the stars (and dust grains) in the 
galaxy differ considerably, such as in a rotating disc. If the dust 
grain velocities are neglected, the mean projected velocity is signif- 
icantly underestimated, and the projected velocity dispersion field 
is completely overestimated. 
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This procedure works equally well when we introduce a cut-off, or 
when a sech or isothermal profile is used instead of the exponential 
profile to describe the vertical distribution of the stars (e.g. Bianchi 
et al. 1996). 

Applying the same technique to generate a random R, we can 
find R by solving the equation 



X = l- 1 + 



R^ 



exp - 



R_ 



(A4) 



Bianchi et al. (1996) invert this transcendental equation numeri- 
cally to find R. However, this equation can be solved exactly by 
means of the Lambert function, yielding 



R = K 



W-i 



X 



(A5) 



The Lambert function, also known as the product log function, is 
generally defined as the inverse of the function w — > /(w) = w e"" , 
and W-i{z) represents the only real branch of this complex func- 
tion besides the principle branch. The Lainbert function can be 
computed in a very efficient way by means of Halley iteration; for 
more information on this function see Corless et al. (1996). 

This paper has been typeset from a TgXy ET^ file prepared by the 
author. 



APPENDIX A: GENERATING POSITIONS FROM AN 
EXPONENTIAL DISC GALAXY 

In this Appendix, we show how the SKIRT code generates a ran- 
dom position from the exponential disc model, i.e. how a random 
position is drawn from the three-dimensional probability density 

p(r)dr oc exp ^^P ^^^^ 

As this probability density is a separable function of R, ip and z, we 
can generate a random position by independently generating each 
of these three coordinates, whereby the generation of an azimuth is 
trivial. To generate a random radius and height, we use the trans- 
formation technique, which says that a random variable x can be 
drawn from a probability density p{x)dx by generating a tmiform 
deviate X and solving the equation 

p{x') dx' / / p{x') dx' (A2) 

-oo / </ —oo 

for X. The calculation of a random z is easy using this principle and 

results in 



z = z» sgn(l + 2X) ln(l - |1 - 2X\). 



(A3) 



